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Abstract 

High resolution, direct numerical simulations of three-dimensional incompressible Navier-Stokes 
equations are carried out to study the energy spectrum in the dissipation range. An energy spectrum 
of the form A(k/kd) Q exp[— (3k/kd\ is confirmed. The possible values of the parameters a and /?, as 
well as their dependence on Reynolds numbers and length scales, are investigated, showing good 
agreement with recent theoretical predictions. A “bottleneck”-type effect is reported at k/kd ~ 4, 
exhibiting a possible transition from near-dissipation to far- dissipation. 
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I. Introduction 


The study of the statistical properties of fluid turbulence at small scales is very difficult from 
both experimental and numerical points of view. In the former, the signal at these scales is exceed- 
ingly weak and consequently it is easily contaminated by the background noise that is naturally 
present in experiments. Resolution of inertial-range and dissipation-range scales at high Reynolds 
numbers by direct numerical simulations (DNS) requires computational capabilities far beyond 
those available today. 

In the last four decades, several functional forms for the energy spectrum in the dissipation 
and far-dissipation ranges have been proposed based on heuristic arguments. The energy spectrum 
in isotropic turbulence can be written as E(k) = e 2 ^ 3 k~ 5 ^ 3 F(k/kd), consistent with Kolmogorov’s 
ideas (Kolmogorov 1941). In the expression above, E(k) = 47rfc 2 (l/2|u(fc)| 2 ), where u(k ) is the 
Fourier transform of the velocity field; kd = e 1 / 4 ^ 3 / 4 is the Kolmogorov dissipation wave number, 
that represents the scale at which the relative importance of the inertial and dissipative terms is 
comparable; c is the mean rate of energy dissipation per unit mass, and v is the kinematic viscosity. 
() denotes an ensemble average. For scales associated with wavenumbers much smaller than kd , the 
energy spectrum becomes essentially independent of the viscosity v. and F(k/kd ) is in that Omit 
the Kolmogorov constant Ck ■ For wavenumbers k > kd-, dissipation is important and the function 
F has a strong dependence on the viscosity through kd . 

A widely proposed form for E(k) in the dissipation range is E(k ) ~ (k/kd) a exp(-f3(k/kd) n ). 
There has been considerable debate as to what the values of the parameters a, / 3 , and n are. 
Townsend (1951) and Novikov (1961) suggested n = 2, based on a Taylor expansion of the velocity 
field at scales much smaller than tj = 1 /kd, the dissipation length scale. This same value is preferred 
by Smith and Reynolds in a more recent work (Smith & Reynolds 1991) in which they compare 
the skewness for the longitudinal velocity derivative obtained from experiments, with the values 
provided by E(k) ~ (k/kd) a ex.p(—0(k/kd) n ) with n = 1,4/3 and 2, and a = -5/3. They conclude 
that n — 2 is the best fit, and that n = 1 is a poor model. Sanada (1992) proposed a spectrum 
with n = 1, but with a sharp change at a transitional wavenumber. His model fit the spectrum 
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from simulations and yielded a skewness of the velocity derivative in agreement with experimental 
values. Manley (1992) argued that the measurements used by Smith and Reynolds did not extend 
much beyond the Kolmogorov scale r/, and accordingly he introduced a cutoff wavenumber for the 
computation of the skewness. With this correction, n = 1 is a better fit for the experimental 
data than n = 2. Kraichnan's direct interaction approximation (DIA) (Kraichnan 1959) predicts 
n = 1 for the far- dissipation range; Orszag (1973) argues n < 1 and Foias et al. conclude n > 1 
(1990) from their analysis. Sreenivasan (1985) found experimentally n = 1 (for 0.5 < k/kd < 1 . 5 ), 
and Chen et al. (1993) observed that the same value would best fit the energy spectrum, in a 
highly resolved DNS at a relatively low Reynolds number. Thus, the value of n = 1 seems to be 
established from theories, experiments, and numerical simulations. In the rest of this report, we 
will concentrate on an energy spectrum of the form 

E(k) ~ (k/k d ) a exp(-/3(k/k d )). (1) 

There are only few studies for determining a and 3. Some phenomenological models (McComb 
1990) suggest a = -5/3, the value that has the attractive property of being consistent with 
Kolmogorov’s scaling for lower wavenumbers. On the other hand, both DIA (Kraichnan 1959) and 
the quasi-normal theories (Orszag 1973) predict a - 3. The recent numerical studies by Kida et 
al. (1992) with data for k/kd < 4 found that a < 0. Kida et al. argue that the negative value 
for a is connected to the fact that in the dissipation range, the energy transfer is dominated by 
nonlocal triads. See also the study by Zhou (1994). Domaradzki (1992) obtains a = -2, deduced 
from a scaling ansatz for energy transfer. On the other hand, Chen et al. (1993) have confirmed 
the prediction of the DIA and the EDQNM for 4 < k/kd < 10- obtaining a = 3.3 and (3 = 7.1. 
Since the Reynolds numbers have been so low, some concerns still remain. Using very large w r ind 
tunnel experiments, recently, Saddoughi & Veeravalli (1994) found 3 = 5.2, for 0.2 < k/kd < 1» 
at R\ ss 600, in accordance with other works (see e.g. Sanada 1992, Kida et al. 1992, Kida & 
Murakami 1990). 
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II. Analysis 


In this note, we concentrate on the determination of a and /?, and their possible dependence 
on the Reynolds number and wavelength, from numerical simulations. A series of direct numerical 
simulations of incompressible Navier-Stokes equations has been carried out in a periodic box with 
sides 27 T. We use a pseudo-spectral code on the CM-5 parallel supercomputer at Los Alamos 
National Laboratory. The resolution ranges from 256 3 (maximum wavenumber k max = 128), to 
512 3 (k maT = 256). The system is initialized with a random velocity field with Gaussian statistics. 
These Fourier velocity amplitudes give an energy spectrum of k 4 exp \-{kjk 0 ) 2 }, with A: 0 ~ 5. If 
the system is then left to evolve, nonlinear interactions will redistribute the energy among all 
available modes. This familiar process delivers energy (in a statistical sense) from larger scales 
to smaller scales, until the excited wavenumbers are high enough for the viscosity to efficiently 
dissipate energy. In order to achieve a steady state, energy is pumped into the system by stirring 
it at the large scales. Specifically, the mode amplitudes for k = 1 and 2 (the largest scales) are 
forced to values consistent with Kolmogorov’s -5/3 law for the energy spectrum. Once statistically 
steady values of R\ and dissipation per unit mass, e, were reached, several snapshots of the energy 
spectrum, spanning several eddy-turnover characteristic times, were averaged and analyzed. 

In Table I, the simulations are sorted by increasing R\ = v 0 \/v, the Taylor microscale Reynolds 
number, where vq is the mean square velocity, and A is the Taylor microscale. The table also displays 
information about how well the dissipation range is resolved (column under “k max / kd' 1 ). Notice 
that although for some of the simulations k max is of the order of kj, for six of the simulations 
k/kd > 5, i.e., fairly deep into the dissipation range. In particular, the third entry in the Table 
(R\ = 38) corresponds to a double-precision, 512 3 DNS, which constitutes to the best knowledge 
of the authors the memory-wise most demanding turbulence simulation to date. 

To investigate the possible values for a and /?, the quantity, din E(k)/dlnk, as a function 
of k/kd- was computed from the energy spectrum, for each of the runs. If the functional form 
suggested in (1) is correct, then din E{k)/d\nk should be close to a straight line. The slope, 0/kd, 
and the intersection of the line with the vertical axis (at din E(k)/ dink = a) were then determined 
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by means of a least-squares fit to restricted ranges of k/k d . 

Figure 1 shows din E(k)/dln k vs k/k d for the run that best resolves the dissipation range 
(run 3 in Table I). Two different slopes are found for the lower and higher wavenumber regions, 
the latter one being steeper. The transition occurs at k/k d ~ 4. These are common features for 
the simulations that penetrate deep into the dissipation range, namely, the two-slope behavior for 
din E(k)/d\nk vs. k/k d , and the transition between the two slopes at k/k d ~ 4. It is possible that 
this shift in slope indicates a change in the nature of the balance of the triadic interactions. We have 
calculated the transfer energy function T{k\p,q) that measures the amount of energy transferred 
to wavenumber k , due to interactions with wavenumbers p and q. A preliminary study of T{k\p.q) 
has not proven to be very revealing. 

A bump at k/k d ~ 4 in Fig. 1 accompanies the occurrence of a change in the energy slope. 
Similar bumps at much smaller k/kd have been ascribed to a “bottleneck” effect (Lohse & Miiller- 
Groeling 1995) in which triadic interactions that can cascade energy out of the region are inhibited 
because the spectrum falls off rapidly as kjk d increases. The energy pileup in our study seems to be 
of a weaker nature than the one reported by Lohse and Miiller-Groeling. In that work, the second- 
order derivative of the energy with respect to k , is non-monotonic. In our study, the third derivative 
is the non-monotonic one. A completely satisfying explanation of this effect is more difficult since 
there is not a clearly conserved quantity associated with the dissipation range, whereas, in the 
inertial range, the energy flux is independent of the wavenumber, i.e., it is conserved. It could 
be conjectured through these observations that the sharp change in the local slope of the energy 
spectrum, will result in an energy pileup close to the transition point, due to the change of the 
nature of the balance of the triadic interactions. 

From the energy spectrum shown in Fig. 1 we computed a and Q as a function of k/k d . The 
result is displayed in Figure 2. Panels a) and b) show a and /?, respectively, as a function of k/kd 
for run 3 in Table I. These values were obtained as the local a and (3 (i.e., local intersection with the 
y-axis and local slope) of din Ejd In k. The values for a and (3 assigned to a particular wavenumber, 
were computed as a least-squares fit to the energy spectrum in the range [k/k d - 0.9 ,k/k d + 0.9]. 
For a we observe a tendency to be negative for lower k and positive for higher wavenumbers. The 
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plot in panel b) shows that 0 is fairly constant at a value between five and six, except at k/k d ~ 4, 
where it drops to about four. For both parameters we observe that a plateau is reached for k/kd < 2 
and k/kd > 8. We note that this plateau is not expected to extend to the inertial range scales since 
that would imply that the energy spectrum would have the same functional form in the dissipation 
and in the inertial ranges. 

Figure 3 shows the dependence of the coefficients a and 0 on R\. Panel a) displays a vs Rx. 
Points indicated by the “x” symbol, correspond to measurements done for k/k d < 3 for all runs in 
Table I; the diamonds correspond to k/k d > 6 for runs 1 through 4. For a, all values obtained from 
the range k/kd < 3 are consistently negative, within the range -1 to -2, in the neighbourhood of 
the Kolmogorov’s value -5/3. The four points obtained in the more dissipative range, are positive, 
in the range 3-6. This two-value behavior of a is a preliminary result. It might be possible that 
the region of positive values for a can be extended beyond R\ ~ 40-50 if these values are evaluated 
from computations that better resolve the dissipation range. In other words, it is possible that the 
negative values for a appear because they are being computed from the near-dissipation range. 

For 0, the values measured in the range k/k d > 6 (diamonds), are consistently higher than 
those measured in the range k/k d < 3. Recently, Sirovich et al. derived the energy spectrum 
(for both inertial range and dissipation range) from the third order structure function (Sirovich. 
Smith & Yakhot 1994). The value of 0 obtained in that work is a function of the skewness and the 
dissipation length. The skewness is believed to be independent of R\ (Kerr 1985), and thus it is 
likely that the behavior of 0 observed in Fig. 2b) is consistent with the derivation in Sirovich et 
al. (1994). 

III. Conclusion 

These values for a and 0 are consistent with those reported by Chen et al. (1993) obtained 
with a DNS with a dissipation range resolution comparable to run number three in Table I, but 
with a lower Reynolds number. The results presented here, and particularly in reference to Fig. 
3, might indicate that the exponential functional form (1) set forth for the energy spectrum in 
the dissipation range, may not be the most appropriate guess for the near-dissipation range, i.e., 
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k/kd < 4. Solidly establishing these parameters as a function of R\ will require next-generation 
computer capabilities. At this point we are faced with the dilemma of having to choose between 
performing DNS with relatively high values for R\ or having a relatively well- resolved dissipation 
range (kmax/kj >> 1) at a relatively low Reynolds number. This study shows a definite trend, and 
it should be taken as an effort to shed some light on these issues by using state-of-the-art DNS. 

Computations were performed on the CM5 at the Advanced Computing Laboratory at Los 
Alamos, and on the SP-1 and SP-2 at the IBM Watson Research Center. Discussions with N. Cao, 
D. Lohse, Z.-S. She, and V. Yakhot are gratefully acknowledged. 
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TABLES 





kmax / kd 

Resolution 

1 

20.9 

9.7 

10 

256 

2 


13.9 

8 

256 

3 


16.6 

12 

512 

4 

39.1 

16.7 

7 

256 

5 

48.4 

19.9 

5.8 

256 

6 

52.9 

23.6 

4.7 

256 

7 

67.9 

33.4 

3.3 

256 

8 



4.4 

512 


TABLE I. Direct numerical simulations on a 27T 3 periodic domain. Resolution ranges are 128 3 


256 3 , and 512 3 . The third entry corresponds to a 512 3 resolution, double precision run. 








FIGURES 



FIG. 1. dlnE/dlnk vs. k/k d for run 3 in Table I. Two slopes are observed in the low and high 
wavenumber range. The energy spectrum E{k ) is computed by averaging the energy spectrum over 
several eddy-turnover times after the system has reached a steady state. 
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FIG. 3. Parameters a a), and Q b), as a function of R \ . These values are calculated by 
performing a least-squares fit to din Ej dink vs. k/kd (see Fig. 1). Points indicated by the “x” 
symbol, correspond to measurements done in the range k/kd < 3 for all runs in Table I; the 
diamonds correspond to k/kd > 6 for runs 1 through 4. 
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